Postmortem spinal cord ALS vs control differential expression
Dr Oliver Ziff
The Francis Crick Institute
Institute of Neurology
University College London, UK
The Francis Crick Institute
Institute of Neurology
University College London, UK
2023-03-01
Fig 4 Postmortem
Volcano
nygc_postmortem_spinal_cord.als.labels <- nygc_postmortem_spinal_cord$res %>% filter(gene_name %in% ALS_genes) %>% top_n(n=5,wt = abs(stat)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.p53.labels <- nygc_postmortem_spinal_cord$res %>% filter(gene_name %in% p53.DDR.repair_gene_set) %>% top_n(n=5,wt = abs(stat)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.labels <- nygc_postmortem_spinal_cord$res %>% filter(!grepl("^AC|^AL", gene_name)) %>% top_n(n=10,wt = abs(stat)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.volcano <- volcano_plot(nygc_postmortem_spinal_cord$res, numerator = "ALS", denomenator = "CTRL",
title = "Postmortem spinal cord", subtitle = "57 Controls vs 214 ALS", ymax = 45, xmax = 2.9, dpi = 300,
gene_labels = c(nygc_postmortem_spinal_cord.als.labels, #nygc_postmortem_spinal_cord.ipsn.up.labels, nygc_postmortem_spinal_cord.ipsn.down.labels,
nygc_postmortem_spinal_cord.labels, nygc_postmortem_spinal_cord.p53.labels))
nygc_postmortem_spinal_cord.volcano
Functional enrichment terms
nygc_postmortem_spinal_cord.gost <- nygc_postmortem_spinal_cord$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>% map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.gprofiles <- nygc_postmortem_spinal_cord.gost$result %>% filter(!str_detect(source, "TF|CORUM|HPA|MIRNA")) %>% arrange( -log10(p_value) ) %>% mutate(query = case_when(query == "query_1" ~ " ALS Down ", query == "query_2" ~ " ALS Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.gprofiles_down <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Down ")
nygc_postmortem_spinal_cord.gprofiles_up <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Up ")
down_list <- c(
"synapse",
"axon",
"neuron differentiation",
"cytoskeleton organization",
"neurogenesis",
"transcription by RNA polymerase II",
"regulation of RNA metabolic process",
"protein modification process")
nygc_postmortem_spinal_cord.gprofiles_filt_down <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"regulation of response to DNA damage stimulus",
"Stabilization of p53",
"programmed cell death",
"endoplasmic reticulum",
"lysosome",
"cellular response to stress",
"response to stress",
"translation",
"protein metabolic process",
"mitochondrion")
nygc_postmortem_spinal_cord.gprofiles_filt_up <- nygc_postmortem_spinal_cord.gprofiles %>% filter(query == " ALS Up ") %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.gprofiles_filt_down, nygc_postmortem_spinal_cord.gprofiles_filt_up) %>% mutate(query = factor(query, levels = c(" ALS Up ", " ALS Down ")))
nygc_postmortem_spinal_cord.gprofiles_filt_combined = nygc_postmortem_spinal_cord.gprofiles_filt_combined %>% mutate(term_name = gsub("regulation of ", "", term_name), term_name = gsub("and","&",term_name), term_name = gsub("cellular ","",term_name), term_name = gsub("response to |transcription by ","",term_name))
nygc_postmortem_spinal_cord.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.gprofiles_filt_combined)# + theme_oz() #theme(legend.position = "none", axis.text = element_text(size=10))
nygc_postmortem_spinal_cord.go_curated
Pathways and transcription factors
net_progeny <- get_progeny(organism = 'human', top = 100)
nygc_postmortem_spinal_cord.als_vs_ctrl.res.matrix <- nygc_postmortem_spinal_cord$res %>% select(gene_name, stat) %>% drop_na(stat) %>% distinct(gene_name, .keep_all = TRUE) %>% tibble_to_matrix(stat, row_names = "gene_name")
nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.contrast_acts <- run_wmean(mat=nygc_postmortem_spinal_cord.als_vs_ctrl.res.matrix, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = source, group2 = source)
nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.pathwayActivity <- ggplot(nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.contrast_acts, aes(x=reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
theme_oz() + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = "Signaling Pathways", x = "", y = "Normalised enrichment ALS vs CTRL") +
stat_pvalue_manual(nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.contrast_acts, label = "p.signif", y.position = 6, xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE)
nygc_postmortem_spinal_cord.als_vs_ctrl_progeny.pathwayActivity
# p53 weights
nygc_postmortem_spinal_cord_progeny.p53weights = net_progeny %>% filter(source == "p53", target %in% nygc_postmortem_spinal_cord$res$gene_name) %>% left_join(select(nygc_postmortem_spinal_cord$res, target=gene_name,stat)) %>%
mutate(color = case_when(weight > 0 & stat > 0 ~ '1', weight > 0 & stat < 0 ~ '2', weight < 0 & stat > 0 ~ '2', weight < 0 & stat < 0 ~ '1', TRUE ~ "3"))
nygc_postmortem_spinal_cord_progeny.p53weights.plot = ggplot(nygc_postmortem_spinal_cord_progeny.p53weights, aes(x = weight, y = stat, color = color)) + geom_point(size = 0.5) +
scale_colour_manual(values = c("red","royalblue3","grey")) +
geom_text_repel(fontface = "italic", data = filter(nygc_postmortem_spinal_cord_progeny.p53weights, (weight > 10 & stat > 2) | (weight > 0 & stat > 6)), aes(label = target), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.1) +
theme_oz() + theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = "p53 pathway weights", x = "p53 pathway gene weight", y = "ALS vs CTRL test statistic") +
geom_vline(xintercept = 0, linetype = 'dotted') + geom_hline(yintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-8,18), ylim = c(-7,10))
nygc_postmortem_spinal_cord_progeny.p53weights.plot
# TFs
net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts <- run_wmean(mat=nygc_postmortem_spinal_cord.als_vs_ctrl.res.matrix, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 10000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>%
mutate(sig = case_when(p_value < 0.06 & abs(score) < 4 ~ "sig", p_value < 0.06 & abs(score) >= 4 ~ "sig_strong", p_value > 0.05 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.up = nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts %>% top_n(score, n = 3) %>% pull(source)
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.down = nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts %>% top_n(-score, n = 3) %>% pull(source)
# Volcano of TFs score vs p_value
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.volcano = ggplot(nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray", "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
labs(y = expression(-log[10]~P~value), x = "Normalised enrichment ALS vs CTRL", title = "Transcription Factors") +
guides(colour = "none") +
theme_oz() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
# scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) + scale_x_continuous(limits = c(-xmax,xmax))
geom_text_repel(fontface = "italic", data = filter(nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts, source %in% c(nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.up, nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.labels.down, "TP53")), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
geom_vline(xintercept = 0, linetype = 'dotted')
nygc_postmortem_spinal_cord.als_vs_ctrl_dorothea.contrast_acts.volcano
Compare postmortem with iPSNs
postmortem_ipsn.als_vs_ctrl = select(nygc_postmortem_spinal_cord$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_als_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.als_vs_ctrl.overlapping = postmortem_ipsn.als_vs_ctrl %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name) # 50
postmortem_ipsn_de_list = list("Postmortem spinal cord: ALS vs CTRL" = nygc_postmortem_spinal_cord$res, "iPSN: ALS vs CTRL" = ipsc_mn_als_datasets.res)
postmortem_ipsn.als_vs_ctrl.scatter = plot_de_correlation(model_list1 = postmortem_ipsn_de_list, model_list2 = postmortem_ipsn_de_list, mutation1 = "iPSN: ALS vs CTRL", mutation2 = "Postmortem spinal cord: ALS vs CTRL", gene_labels = c(postmortem_ipsn.als_vs_ctrl.overlapping, "CFAP410"), plot_p_value = TRUE)
postmortem_ipsn.als_vs_ctrl.scatter
Compare mutations in postmortem
# heatmap of correlation coefficients
postmortem_mutations_stat = select(nygc_postmortem_spinal_cord.sporadic$res, gene_id, Sporadic = stat) %>%
left_join(select(nygc_postmortem_spinal_cord.c9orf72$res, gene_id, C9orf72 = stat)) %>%
left_join(select(nygc_postmortem_spinal_cord.sod1$res, gene_id, SOD1 = stat)) %>%
left_join(select(nygc_postmortem_spinal_cord.fus$res, gene_id, FUS = stat)) %>%
drop_na() %>%
column_to_rownames("gene_id")
cormat <- round(cor(postmortem_mutations_stat),2)
melted_cormat <- melt(cormat)
get_lower_tri<-function(cormat){ # Get lower triangle of the correlation matrix
cormat[upper.tri(cormat)] <- NA
return(cormat)
}
get_upper_tri <- function(cormat){ # Get upper triangle of the correlation matrix
cormat[lower.tri(cormat)]<- NA
return(cormat)
}
upper_tri <- get_upper_tri(cormat)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cormat){ # Use correlation between variables as distance
dd <- as.dist((1-cormat)/2)
hc <- hclust(dd)
cormat <-cormat[hc$order, hc$order]
}
cormat <- reorder_cormat(cormat)
upper_tri <- get_upper_tri(cormat)
melted_cormat <- melt(upper_tri, na.rm = TRUE) # Melt the correlation matrix
postmortem.dge_correlation_ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white") + scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, limit = c(-1,1), space = "Lab", name="Pearson\nCorrelation") +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 3) +
theme_minimal() + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), panel.grid.major = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.ticks = element_blank(),
legend.position = "top", legend.direction = "horizontal", legend.title = element_text(size = 8), axis.text=element_text(size=7.5)) + # coord_fixed()+
guides(fill = guide_colorbar(barwidth = 7, barheight = 1, title.position = "top", title.hjust = 0.5))
postmortem.dge_correlation_ggheatmap
p53 pathway and TP53 TF
nygc_postmortem_spinal_cord_mn_mutations_list = list("Sporadic" = select(nygc_postmortem_spinal_cord.sporadic$res,gene_name,stat), "C9orf72" = select(nygc_postmortem_spinal_cord.c9orf72$res,gene_name,stat), "SOD1" = select(nygc_postmortem_spinal_cord.sod1$res,gene_name,stat), "FUS" = select(nygc_postmortem_spinal_cord.fus$res,gene_name,stat))
mutations <- names(nygc_postmortem_spinal_cord_mn_mutations_list)
nygc_postmortem_spinal_cord_mn_mutations_list = map(nygc_postmortem_spinal_cord_mn_mutations_list, ~{drop_na(.x)})
nygc_postmortem_spinal_cord_mn_mutations_list = map(nygc_postmortem_spinal_cord_mn_mutations_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
nygc_postmortem_spinal_cord_mn_mutations_list = map(nygc_postmortem_spinal_cord_mn_mutations_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})
# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
nygc_postmortem_spinal_cord_mn_mutations_list.progeny_scores = map_df(nygc_postmortem_spinal_cord_mn_mutations_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean')
nygc_postmortem_spinal_cord_mn_mutations_list.p53.progeny_scores = nygc_postmortem_spinal_cord_mn_mutations_list.progeny_scores %>% filter(source == "p53") %>% mutate(decoupleR = "p53 pathway")
# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
nygc_postmortem_spinal_cord_mn_mutations_list.dorothea_scores = map_df(nygc_postmortem_spinal_cord_mn_mutations_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean')
nygc_postmortem_spinal_cord_mn_mutations_list.TP53.dorothea_scores = nygc_postmortem_spinal_cord_mn_mutations_list.dorothea_scores %>% filter(source == "TP53") %>% mutate(decoupleR = "TP53 TF")
nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53 = bind_rows(nygc_postmortem_spinal_cord_mn_mutations_list.p53.progeny_scores, nygc_postmortem_spinal_cord_mn_mutations_list.TP53.dorothea_scores) %>%
mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)
nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53.plot <- ggplot(nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53, aes(x=reorder(mutation, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) + #scale_fill_npg() +
theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 30, hjust = 1), plot.title = element_text(hjust = 0.5), legend.position = "none") +
geom_hline(yintercept = 0, linetype = 3) +
labs(x = "", y = "Normalised enrichment") + #ylim(c(-12,12)) +
stat_pvalue_manual(nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53, label = "p.signif", y.position = 6, xmin = "mutation", xmax = NULL, size = 4, hide.ns = TRUE) +
facet_wrap(~decoupleR, scales = "free")
nygc_postmortem_spinal_cord_mn_mutations_list.p53.TP53.plot
Table S7 postmortem & ipsn overlap DGE
Overlapping differentially expressed genes in ALS vs CTRL in both postmortem spinal cord and iPSNs
postmortem_ipsc_overlap.res =
mutate(select(nygc_postmortem_spinal_cord$res, gene_name, gene_id, padj.postmortem = padj, stat.postmortem = stat, lfc.postmortem = log2FoldChange), direction.postmortem = case_when(padj.postmortem < 0.05 & lfc.postmortem > 0 ~ "up", padj.postmortem < 0.05 & lfc.postmortem < 0 ~ "down", TRUE ~ "non-significant")) %>%
full_join(mutate(select(ipsc_mn_als_datasets.res, gene_name, gene_id, padj.ipsn = padj, stat.ipsn = stat, lfc.ipsn = log2FoldChange), direction.ipsn = case_when(padj.ipsn < 0.05 & lfc.ipsn > 0 ~ "up", padj.ipsn < 0.05 & lfc.ipsn < 0 ~ "down", TRUE ~ "non-significant"))) %>%
select(gene_name, gene_id, everything()) %>%
filter(padj.ipsn < 0.05, padj.postmortem < 0.05, ((stat.ipsn > 0 & stat.postmortem > 0) | (stat.ipsn < 0 & stat.postmortem < 0)))
kbl(postmortem_ipsc_overlap.res) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), fixed_thead = T) %>% scroll_box(width = "100%", height = "500px")| gene_name | gene_id | padj.postmortem | stat.postmortem | lfc.postmortem | direction.postmortem | padj.ipsn | stat.ipsn | lfc.ipsn | direction.ipsn |
|---|---|---|---|---|---|---|---|---|---|
| ISCU | ENSG00000136003 | 0.0000000 | 8.311838 | 0.4136826 | up | 0.0118787 | 4.469863 | 0.1030001 | up |
| RRM2B | ENSG00000048392 | 0.0000000 | 8.233753 | 0.7036306 | up | 0.0167379 | 4.352172 | 0.2020891 | up |
| PAQR4 | ENSG00000162073 | 0.0000000 | -7.967077 | -0.5929364 | down | 0.0397368 | -4.022565 | -0.1423603 | down |
| RNASEL | ENSG00000135828 | 0.0000000 | 6.810334 | 0.3033616 | up | 0.0032664 | 5.089660 | 0.5962462 | up |
| PTCHD4 | ENSG00000244694 | 0.0000000 | 6.182531 | 0.5597501 | up | 0.0086407 | 4.663660 | 0.3292627 | up |
| PLIN5 | ENSG00000214456 | 0.0000001 | -5.759645 | -0.7603759 | down | 0.0328980 | -4.107068 | -0.4650748 | down |
| FBXO22 | ENSG00000167196 | 0.0000067 | 4.891648 | 0.2366739 | up | 0.0461030 | 3.967703 | 0.1144354 | up |
| SOD2-OT1 | ENSG00000285427 | 0.0000255 | -4.587429 | -0.4386346 | down | 0.0170913 | -4.337828 | -0.2587170 | down |
| SESN1 | ENSG00000080546 | 0.0001084 | 4.240309 | 0.2871041 | up | 0.0163788 | 4.377748 | 0.1711968 | up |
| MTCO2P12 | ENSG00000229344 | 0.0001459 | -4.165748 | -0.7582162 | down | 0.0092457 | -4.621477 | -1.4451766 | down |
| OLIG1 | ENSG00000184221 | 0.0009450 | -3.663960 | -0.3287542 | down | 0.0328980 | -4.101139 | -0.7260811 | down |
| OLIG2 | ENSG00000205927 | 0.0010312 | -3.638942 | -0.3100081 | down | 0.0297222 | -4.154218 | -0.7014318 | down |
| CDKN1A | ENSG00000124762 | 0.0038023 | 3.243198 | 0.5623347 | up | 0.0198326 | 4.283226 | 0.3542375 | up |
| CROCC2 | ENSG00000226321 | 0.0088169 | -2.961274 | -1.0379353 | down | 0.0253547 | -4.206046 | -1.5582071 | down |
| RPL29P11 | ENSG00000224858 | 0.0115688 | -2.864505 | -0.6140480 | down | 0.0397368 | -4.020057 | -0.9411368 | down |
| AC025171.1 | ENSG00000177738 | 0.0283717 | -2.524708 | -0.1924196 | down | 0.0077084 | -4.755944 | -0.3937060 | down |
| MTCO1P12 | ENSG00000237973 | 0.0356193 | -2.432034 | -0.5159713 | down | 0.0000000 | -8.158533 | -1.6010349 | down |
Scatter postmortem & iPSN subgroups
postmortem_ipsc_mn_mutations_list = list("Sporadic iPSN" = ipsc_mn_sporadic_datasets.res, "C9orf72 iPSN" = ipsc_mn_c9orf72_datasets.res, "SOD1 iPSN" = ipsc_mn_sod1_datasets.res, "FUS iPSN" = ipsc_mn_fus_datasets.res, "Sporadic postmortem" = nygc_postmortem_spinal_cord$res, "C9orf72 postmortem" = nygc_postmortem_spinal_cord.c9orf72$res, "SOD1 postmortem" = nygc_postmortem_spinal_cord.sod1$res, "FUS postmortem" = nygc_postmortem_spinal_cord.fus$res)
postmortem_ipsn.sporadic.join = select(nygc_postmortem_spinal_cord$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_sporadic_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.sporadic.overlapping = postmortem_ipsn.sporadic.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsn.c9orf72.join = select(nygc_postmortem_spinal_cord.c9orf72$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_c9orf72_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.c9orf72.overlapping = postmortem_ipsn.c9orf72.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsn.sod1.join = select(nygc_postmortem_spinal_cord.sod1$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_sod1_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.sod1.overlapping = postmortem_ipsn.sod1.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsn.fus.join = select(nygc_postmortem_spinal_cord.fus$res, gene_name, gene_id, nygc.stat = stat, nygc.padj = padj) %>% left_join(select(ipsc_mn_fus_datasets.res, gene_name, gene_id, ipsn.stat = stat, ipsn.padj = padj))
postmortem_ipsn.fus.overlapping = postmortem_ipsn.fus.join %>% filter(nygc.padj < 0.05, ipsn.padj < 0.05, (nygc.stat > 0 & ipsn.stat > 0) | (nygc.stat < 0 & ipsn.stat < 0)) %>% top_n(n = 20, wt = abs(nygc.stat)) %>% pull(gene_name)
postmortem_ipsc_mn_mutations_list_t_stat_plot <- plot_grid(
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "Sporadic iPSN", mutation2 = "Sporadic postmortem", gene_labels = postmortem_ipsn.sporadic.overlapping, plot_p_value = TRUE),
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "C9orf72 iPSN", mutation2 = "C9orf72 postmortem", gene_labels = c(postmortem_ipsn.c9orf72.overlapping, "C9orf72"), plot_p_value = TRUE),
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "SOD1 iPSN", mutation2 = "SOD1 postmortem", gene_labels = c(postmortem_ipsn.sod1.overlapping, "SOD1"), plot_p_value = TRUE),
plot_de_correlation(model_list1 = postmortem_ipsc_mn_mutations_list, model_list2 = postmortem_ipsc_mn_mutations_list, mutation1 = "FUS iPSN", mutation2 = "FUS postmortem", gene_labels = c(postmortem_ipsn.fus.overlapping, "FUS"), plot_p_value = TRUE),
ncol = 2, nrow = 2)# & xlim(-8,8) & ylim(-8,8)
postmortem_ipsc_mn_mutations_list_t_stat_plot
Fig S14 Compare mutations postmortem
Compare sporadic, C9orf72, SOD1, FUS
Volcano
postmortem_mutants_bind_datasets.res = bind_rows(mutate(nygc_postmortem_spinal_cord.sporadic$res, mutation = "Sporadic"), mutate(nygc_postmortem_spinal_cord.c9orf72$res, mutation = "C9orf72"), mutate(nygc_postmortem_spinal_cord.sod1$res, mutation = "SOD1"), mutate(nygc_postmortem_spinal_cord.fus$res, mutation = "FUS"))
postmortem_mutants_join_datasets.res = full_join(select(nygc_postmortem_spinal_cord.c9orf72$res, gene_name, padj.c9orf72 = padj, log2FoldChange.c9orf72 = log2FoldChange),
select(nygc_postmortem_spinal_cord.sporadic$res, gene_name, padj.sals = padj, log2FoldChange.sals = log2FoldChange)) %>%
full_join(select(nygc_postmortem_spinal_cord.fus$res, gene_name, padj.fus = padj, log2FoldChange.fus = log2FoldChange)) %>%
full_join(select(nygc_postmortem_spinal_cord.sod1$res, gene_name, padj.sod1 = padj, log2FoldChange.sod1 = log2FoldChange)) #%>% full_join(select(postmortem_vcp_datasets$res, gene_name, padj.vcp = padj, log2FoldChange.vcp = log2FoldChange))
postmortem_mutants_join_datasets.res.labels = postmortem_mutants_join_datasets.res %>% filter(c(padj.c9orf72 < 0.05 & padj.fus < 0.05 & padj.sod1 < 0.05) |
c(padj.fus < 0.05 & padj.sod1 < 0.05) |
c(padj.fus < 0.05 & padj.c9orf72 < 0.05)) %>% pull(gene_name)
nygc_postmortem_spinal_cord.fus.labels <- nygc_postmortem_spinal_cord.fus$res %>% filter(-log10(padj)>3, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.fus.volcano <- volcano_plot(nygc_postmortem_spinal_cord.fus$res, "FUS", numerator = "FUS", denomenator = "CTRL", subtitle = "57 Controls vs 2 FUS", ymax = 15, xmax = 4, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.fus.labels))
nygc_postmortem_spinal_cord.fus.volcano
nygc_postmortem_spinal_cord.c9orf72.labels <- nygc_postmortem_spinal_cord.c9orf72$res %>% filter(-log10(padj)>10, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.c9orf72.volcano <- volcano_plot(nygc_postmortem_spinal_cord.c9orf72$res, "C9orf72", numerator = "C9orf72", denomenator = "CTRL", subtitle = "57 Controls vs 36 C9orf72", ymax = 30, xmax = 3, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.c9orf72.labels))
nygc_postmortem_spinal_cord.c9orf72.volcano
nygc_postmortem_spinal_cord.sod1.labels <- nygc_postmortem_spinal_cord.sod1$res %>% filter(-log10(padj)>5, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.sod1.volcano <- volcano_plot(nygc_postmortem_spinal_cord.sod1$res, "SOD1", numerator = "SOD1", denomenator = "CTRL", subtitle = "57 Controls vs 5 SOD1", ymax = 17, xmax = 3, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.sod1.labels))
nygc_postmortem_spinal_cord.sod1.volcano
nygc_postmortem_spinal_cord.sporadic.labels <- nygc_postmortem_spinal_cord.sporadic$res %>% filter(-log10(padj)>15, gene_name %in% c(p53.DDR.repair_gene_set, ALS_genes)) %>% distinct(gene_name) %>% pull(gene_name)
nygc_postmortem_spinal_cord.sporadic.volcano <- volcano_plot(nygc_postmortem_spinal_cord.sporadic$res, "Sporadic", numerator = "Sporadic", denomenator = "CTRL", subtitle = "57 Controls vs 161 Sporadic", ymax = 35, xmax = 3, dpi = 300, gene_labels = c(nygc_postmortem_spinal_cord.sporadic.labels))
nygc_postmortem_spinal_cord.sporadic.volcano
Functional enrichment terms
# c9orf72
nygc_postmortem_spinal_cord.c9orf72.dge_genes <- nygc_postmortem_spinal_cord.c9orf72$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>% map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.c9orf72.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|MIRNA|HP")) %>% arrange(-log10(p_value)) %>%
mutate(query = case_when(query == "query_1" ~ " C9orf72 Down ", query == "query_2" ~ " C9orf72 Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_down <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Down ")
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_up <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Up ")
down_list <- c(
"regulation of RNA metabolic process",
"ribonucleotide binding",
"cellular localization",
"glutamatergic synapse",
"cytoskeleton organization",
"protein modification process",
"synaptic signaling",
"axon",
"nervous system development",
"synapse",
"protein binding")
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"intrinsic apoptotic signaling pathway by p53 class mediator",
"cellular response to DNA damage stimulus",
"Nonsense-Mediated Decay (NMD)",
"apoptotic process",
"programmed cell death",
"cell death",
"mitochondrion",
"cytoplasmic translation",
"response to stress",
"extracellular exosome",
"protein metabolic process",
"protein binding")
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt %>% filter(query == " C9orf72 Up ") %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_up) %>%
mutate(query = factor(query, levels = c(" C9orf72 Up ", " C9orf72 Down ")))
nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined$term_name <- gsub("intrinsic apoptotic signaling pathway by p53 class mediator", "apoptotic signaling by p53", nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined$term_name) %>% gsub("regulation of ", "", .) %>% gsub("response to ","",.) %>% gsub(" \\(NMD\\)","",.)
nygc_postmortem_spinal_cord.c9orf72.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.c9orf72.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.c9orf72.go_curated
# FUS
nygc_postmortem_spinal_cord.fus.dge_genes <- nygc_postmortem_spinal_cord.fus$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>% map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.fus.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|HP|MIRNA")) %>% arrange(-log10(p_value)) %>%
mutate(query = case_when(query == "query_1" ~ " FUS Down ", query == "query_2" ~ " FUS Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.fus.dge_gprofiles_down <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Down ")
nygc_postmortem_spinal_cord.fus.dge_gprofiles_up <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Up ")
down_list <- c(
"oligodendrocyte differentiation",
"cellular localization",
"protein modification process",
"cytoskeleton organization",
"localization",
"axon",
"cytoskeleton",
"protein binding",
"neuron development",
"synapse",
"nervous system development")
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Down ") %>% filter(term_name %in% down_list)
up_list <- c(
# "mitochondrion",
"response to heat",
# "regulation of programmed cell death",
# "response to stress",
"programmed cell death",
"apoptotic process",
"cell death",
"Protein processing in endoplasmic reticulum",
"cellular response to stress",
"response to unfolded protein",
"cytoplasmic translation",
"protein metabolic process",
"protein folding",
# "protein binding",
"extracellular exosome")
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt %>% filter(query == " FUS Up ") %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_up) %>%
mutate(query = factor(query, levels = c(" FUS Up ", " FUS Down ")))
nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined$term_name <- gsub("Protein processing in e", "E", nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined$term_name) %>% gsub("positive regulation of ","",.) %>% gsub("regulation of ","",.)
nygc_postmortem_spinal_cord.fus.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.fus.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.fus.go_curated
# sod1
nygc_postmortem_spinal_cord.sod1.dge_genes <- nygc_postmortem_spinal_cord.sod1$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>% map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.sod1.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|HP|MIRNA")) %>% arrange(-log10(p_value)) %>%
mutate(query = case_when(query == "query_1" ~ " SOD1 Down ", query == "query_2" ~ " SOD1 Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_down <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Down ")
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_up <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Up ")
down_list <- c(
"oligodendrocyte differentiation",
# "cellular localization",
"protein modification process",
"cytoskeleton organization",
# "localization",
"axon",
"cytoskeleton",
"protein binding",
"neuron development",
"synapse",
"nervous system development")
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"cellular response to DNA damage stimulus",
# "p53-Independent DNA Damage Response",
"regulation of immune response",
"cytokine production",
"programmed cell death",
"cell death",
"translation",
"endoplasmic reticulum",
"response to stress",
"protein localization",
# "mitochondrial membrane",
"cellular localization",
"lysosome",
"mitochondrion",
# "extracellular exosome",
"protein metabolic process")
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt %>% filter(query == " SOD1 Up ") %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_up) %>%
mutate(query = factor(query, levels = c(" SOD1 Up ", " SOD1 Down ")))
nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined$term_name <- gsub("Protein processing in e", "E", nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined$term_name) %>% gsub("positive regulation of ","",.) %>% gsub("regulation of ","",.) %>% gsub("cellular response to ","",.)
nygc_postmortem_spinal_cord.sod1.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.sod1.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.sod1.go_curated
# sporadic
nygc_postmortem_spinal_cord.sporadic.dge_genes <- nygc_postmortem_spinal_cord.sporadic$res %>% filter(padj < 0.05) %>% group_split(log2FoldChange > 0) %>% map("gene_id") %>% gost()
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt <- nygc_postmortem_spinal_cord.sporadic.dge_genes$result %>% filter(!str_detect(source, "TF|CORUM|HPA|MIRNA")) %>% arrange(-log10(p_value)) %>%
mutate(query = case_when(query == "query_1" ~ " Sporadic Down ", query == "query_2" ~ " Sporadic Up "), term_name = as.factor(term_name))
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_down <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Down ")
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_up <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Up ")
down_list <- c(
"histone modification",
"cytoskeleton organization",
"neurogenesis",
"nervous system development",
"organelle organization",
"positive regulation of cellular process",
"protein modification process",
"DNA-templated transcription",
"protein binding")
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_down <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Down ") %>% filter(term_name %in% down_list)
up_list <- c(
"cellular response to DNA damage stimulus",
"Stabilization of p53",
# "p53-Dependent G1 DNA Damage Response",
"G1/S DNA Damage Checkpoints",
"Metabolism of RNA",
"regulation of immune system process",
"cell death",
"protein localization",
"endoplasmic reticulum",
"lysosome",
"response to stress",
"cytoplasmic translation",
"protein metabolic process",
"mitochondrion",
"protein binding")
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_up <- nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt %>% filter(query == " Sporadic Up ") %>% filter(term_name %in% up_list)
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined <- bind_rows(nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_down, nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_up) %>%
mutate(query = factor(query, levels = c(" Sporadic Up ", " Sporadic Down ")))
nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined$term_name <- gsub("Protein processing in e", "E", nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined$term_name) %>% gsub("positive regulation of ","",.) %>% gsub("regulation of ","",.) %>% gsub("cellular response to ","",.)
nygc_postmortem_spinal_cord.sporadic.go_curated <- curated_profiles(gprofiles = nygc_postmortem_spinal_cord.sporadic.dge_gprofiles_filt_combined) + theme(legend.position = "none")
nygc_postmortem_spinal_cord.sporadic.go_curated
Upset plot
### Upset plot https://github.com/const-ae/ggupset
postmortem_mutants_bind_datasets.res = bind_rows(mutate(nygc_postmortem_spinal_cord.sporadic$res, mutation = "Sporadic"), mutate(nygc_postmortem_spinal_cord.c9orf72$res, mutation = "C9orf72"),mutate(nygc_postmortem_spinal_cord.sod1$res, mutation = "SOD1"), mutate(nygc_postmortem_spinal_cord.fus$res, mutation = "FUS"))
postmortem_mutants_bind_datasets.res.up.upset = upset_plot(filter(postmortem_mutants_bind_datasets.res, stat > 0), overlap_by = "mutation", overlap_level = "gene_id", split_direction = FALSE, ymax = 2800, order = "freq", order_group_list = c("Sporadic", "C9orf72", "SOD1", "FUS"), title = "Upregulated", label_colours = c("dodgerblue2", "firebrick2", "forestgreen", "gold2"))
postmortem_mutants_bind_datasets.res.up.upset
postmortem_mutants_bind_datasets.res.down.upset = upset_plot(filter(postmortem_mutants_bind_datasets.res, stat < 0), overlap_by = "mutation", overlap_level = "gene_id", split_direction = FALSE, ymax = 2500, order = "freq", order_group_list = c("Sporadic", "C9orf72", "SOD1", "FUS"), title = "Downregulated", label_colours = c("dodgerblue2", "firebrick2", "forestgreen", "gold2"))
postmortem_mutants_bind_datasets.res.down.upset
Pathways
postmortem_mutations_list = list("Sporadic" = select(nygc_postmortem_spinal_cord.sporadic$res,gene_name,stat), "C9orf72" = select(nygc_postmortem_spinal_cord.c9orf72$res,gene_name,stat), "SOD1" = select(nygc_postmortem_spinal_cord.sod1$res,gene_name,stat), "FUS" = select(nygc_postmortem_spinal_cord.fus$res,gene_name,stat))
mutations <- names(postmortem_mutations_list)
postmortem_mutations_list = map(postmortem_mutations_list, ~{drop_na(.x)})
postmortem_mutations_list = map(postmortem_mutations_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
postmortem_mutations_list = map(postmortem_mutations_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})
# net_progeny <- get_progeny(organism = 'human', top = 100)
postmortem_mutations_list.progeny_scores = map_df(postmortem_mutations_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation, source = as.factor(source), source = fct_reorder(source, score), mutation = factor(mutation, levels = c("Sporadic","C9orf72","FUS","SOD1")))
postmortem_mutations_list.progeny_scores.plot <- ggplot(postmortem_mutations_list.progeny_scores, aes(x=mutation, y = score, fill = mutation)) +
geom_bar(stat='identity', position = "dodge2") +
scale_fill_npg() + # scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
labs(x = "", y = "Signalling Pathways\nALS vs CTRL\nNormalised enrichment") +
theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90, size = 6), plot.title = element_text(hjust = 0.5), legend.title = element_blank(), legend.position = "none") +
geom_hline(yintercept = 0, linetype = 3) +
facet_grid(~source) +
stat_pvalue_manual(mutate(postmortem_mutations_list.progeny_scores, y.position = case_when(mutation=="Sporadic"~6.5,mutation=="C9orf72"~7,TRUE~7)), label = "p.signif", xmin = "mutation", xmax = NULL, size = 4, hide.ns = TRUE)
postmortem_mutations_list.progeny_scores.plot
Transcription factors
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
postmortem_mutations_list.dorothea_scores = map_df(postmortem_mutations_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation) %>%
mutate(sig = case_when(p_value < 0.1 & abs(score) < 4 ~ "sig", p_value < 0.1 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.1 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction), mutation = factor(mutation, levels = c("Sporadic","C9orf72","SOD1","FUS")))
# Volcano of TFs score vs p_value
postmortem_mutations_merged_dorothea.contrast_acts.volcano = ggplot(postmortem_mutations_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray", "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
labs(y = expression(-log[10]~P~value), x = "Transcription Factors ALS vs CTRL Normalised enrichment") +
guides(colour = "none") +
theme_oz() + #theme(plot.title = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "C9orf72"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "FUS"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "SOD1"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
geom_text_repel(fontface = "italic", data = top_n(filter(postmortem_mutations_list.dorothea_scores, p_value < 0.05, mutation == "Sporadic"),n = 3, wt = abs(score)), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.3) +
geom_text_repel(fontface = "italic", data = filter(postmortem_mutations_list.dorothea_scores,source %in% c("TP53")), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 2.8) +
geom_vline(xintercept = 0, linetype = 'dotted') +
facet_wrap(~mutation, scales = "free", ncol = 4) & ylim(0,2.8)
postmortem_mutations_merged_dorothea.contrast_acts.volcano
Fig S15 TDP43 vs non TDP43 proteinopathy
iPSMN
load(file = here(proj_path, "expression/deseq2/ipsc_mn_tdp43_status.RData"))
ipsc_mn_tdp43_pos_neg_list = list("TDP-43 ALS" = select(ipsc_mn_tdp43_positive_datasets$res,gene_name,stat), "non-TDP-43 ALS" = select(ipsc_mn_tdp43_negative_datasets$res,gene_name,stat))
ipsc_mn_tdp43_pos_neg_list = map(ipsc_mn_tdp43_pos_neg_list, ~{drop_na(.x)})
ipsc_mn_tdp43_pos_neg_list = map(ipsc_mn_tdp43_pos_neg_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
ipsc_mn_tdp43_pos_neg_list = map(ipsc_mn_tdp43_pos_neg_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})
# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
ipsc_mn_tdp43_pos_neg_list.progeny_scores = map_df(ipsc_mn_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)
ipsc_mn_tdp43_positive_datasets_progeny.pathwayActivity <- ggplot(ipsc_mn_tdp43_pos_neg_list.progeny_scores, aes(x=reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = "iPSMN signaling", x = "", y = "ALS vs CTRL\nNormalised enrichment") +
stat_pvalue_manual(mutate(ipsc_mn_tdp43_pos_neg_list.progeny_scores, y.position = case_when(mutation == "non-TDP-43 ALS" & source %in% c("PI3K","TGFb")~ 4, mutation == "TDP-43 ALS" & source %in% c("MAPK") ~ 14, mutation == "TDP-43 ALS"~15, mutation == "non-TDP-43 ALS"~4.5)), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) +
facet_wrap(~mutation, scales = "free")
ipsc_mn_tdp43_positive_datasets_progeny.pathwayActivity
# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
ipsc_mn_tdp43_pos_neg_list.dorothea_scores = map_df(ipsc_mn_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))
ipsc_mn_tdp43_positive_datasets_dorothea.contrast_acts.volcano = ggplot(ipsc_mn_tdp43_pos_neg_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray", "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
labs(y = expression(-log[10]~P~value), x = "Normalised enrichment ALS vs CTRL", title = "iPSMN TFs") +
guides(colour = "none") +
theme_oz() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
geom_text_repel(fontface = "italic", data = filter(ipsc_mn_tdp43_pos_neg_list.dorothea_scores, source %in% "TP53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +
geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8)) +
facet_wrap(~mutation)
ipsc_mn_tdp43_positive_datasets_dorothea.contrast_acts.volcano
Postmortem
load(file = here(proj_path, "expression/deseq2/nygc_postmortem_spinal_cord_tdp43_status.RData"))
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = list("TDP-43 ALS" = select(nygc_postmortem_spinal_cord_tdp43_positive_datasets$res,gene_name,stat), "non-TDP-43 ALS" = select(nygc_postmortem_spinal_cord_tdp43_negative_datasets$res,gene_name,stat))
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = map(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{drop_na(.x)})
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = map(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
nygc_postmortem_spinal_cord_tdp43_pos_neg_list = map(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})
# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
nygc_postmortem_spinal_cord_tdp43_pos_neg_list.progeny_scores = map_df(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)
nygc_postmortem_spinal_cord_tdp43_positive_datasets_progeny.pathwayActivity <- ggplot(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.progeny_scores, aes(x=reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = "Postmortem signaling", x = "", y = "ALS vs CTRL\nNormalised enrichment") +
stat_pvalue_manual(mutate(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.progeny_scores, y.position = case_when(mutation == "non-TDP-43 ALS" & source == "TNFa" ~ 6, mutation == "TDP-43 ALS" & source %in% c("TNFa","NFkB") ~ 5, mutation == "TDP-43 ALS"~5.5, mutation == "non-TDP-43 ALS"~6.5)), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) +
facet_wrap(~mutation, scales = "free")
nygc_postmortem_spinal_cord_tdp43_positive_datasets_progeny.pathwayActivity
# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
nygc_postmortem_spinal_cord_tdp43_pos_neg_list.dorothea_scores = map_df(nygc_postmortem_spinal_cord_tdp43_pos_neg_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))
nygc_postmortem_spinal_cord_tdp43_positive_datasets_dorothea.contrast_acts.volcano = ggplot(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray", "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
labs(y = expression(-log[10]~P~value), x = "Normalised enrichment ALS vs CTRL", title = "Postmortem TFs") +
guides(colour = "none") +
theme_oz() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
geom_text_repel(fontface = "italic", data = filter(nygc_postmortem_spinal_cord_tdp43_pos_neg_list.dorothea_scores, source %in% "TP53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +
geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8)) +
facet_wrap(~mutation, scales = "free")
nygc_postmortem_spinal_cord_tdp43_positive_datasets_dorothea.contrast_acts.volcano
TDP43 knockdown
Liu TDP43+ vs - postmortem TDP43 knockdown models
neuronal_nuclei_tdp43_liu = readRDS(here(proj_path,"expression/deseq2/neuronal_nuclei_tdp43_liu.rds"))
tdp43_kd_datasets = readRDS(here(proj_path,"expression/deseq2/tdp43_kd_datasets.rds"))
ipsc_mn_tdp43_neuron_kd_list = list("Neuronal Nuclei Brain" = select(neuronal_nuclei_tdp43_liu$res,gene_name,stat),
"TDP-43 knockdown cells" = select(tdp43_kd_datasets$res,gene_name,stat))
ipsc_mn_tdp43_neuron_kd_list = map(ipsc_mn_tdp43_neuron_kd_list, ~{drop_na(.x)})
ipsc_mn_tdp43_neuron_kd_list = map(ipsc_mn_tdp43_neuron_kd_list, ~{distinct(.x, gene_name, .keep_all = TRUE)})
ipsc_mn_tdp43_neuron_kd_list = map(ipsc_mn_tdp43_neuron_kd_list, ~{tibble_to_matrix(.x, stat, row_names = "gene_name")})
# PROGENy
# net_progeny <- get_progeny(organism = 'human', top = 100)
ipsc_mn_tdp43_neuron_kd_list.progeny_scores = map_df(ipsc_mn_tdp43_neuron_kd_list, ~{run_wmean(.x, net=net_progeny, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""), group1 = mutation, group2 = mutation)
ipsc_mn_tdp43_neuron_kd_datasets_progeny.pathwayActivity <- ggplot(ipsc_mn_tdp43_neuron_kd_list.progeny_scores, aes(x=reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = "TDP-43 depletion", x = "", y = "TDP-43 depleted vs TDP-43 retained\nNormalised enrichment") +
stat_pvalue_manual(mutate(ipsc_mn_tdp43_neuron_kd_list.progeny_scores, y.position = case_when(mutation == "Neuronal Nuclei Brain" & source == "p53" ~ 3.5, mutation == "TDP-43 knockdown cells" & source %in% c("TNFa","VEGF") ~ 15, mutation == "TDP-43 knockdown cells"~16, mutation == "Neuronal Nuclei Brain"~4)), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE) +
facet_wrap(~mutation, scales = "free")
ipsc_mn_tdp43_neuron_kd_datasets_progeny.pathwayActivity
# Dorothea
# net_dorothea <- get_dorothea(organism='human', levels=c('A', 'B', 'C'))
ipsc_mn_tdp43_neuron_kd_list.dorothea_scores = map_df(ipsc_mn_tdp43_neuron_kd_list, ~{run_wmean(.x, net=net_dorothea, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5)}, .id = "mutation") %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))
ipsc_mn_tdp43_neuron_kd_datasets_dorothea.contrast_acts.volcano = ggplot(ipsc_mn_tdp43_neuron_kd_list.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray", "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
labs(y = expression(-log[10]~P~value), x = "Normalised enrichment TDP-43 depleted vs TDP-43 positive", title = "TDP-43 depletion TFs") +
guides(colour = "none") +
theme_oz() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
geom_text_repel(fontface = "italic", data = filter(ipsc_mn_tdp43_neuron_kd_list.dorothea_scores, source %in% "TP53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +
geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8)) +
facet_wrap(~mutation, scales = "free")
# ipsc_mn_tdp43_neuron_kd_datasets_dorothea.contrast_acts.volcanoTDP43 overexpression
Mouse model from Maor-nof et al
maor_nof_cell_tdp_vs_ctrl = read_csv(here(proj_path, "expression/maor_nof_cell_tdp/GSE162048_D1CON_vs_D1TDP-ExpDiff.csv")) %>% rename(gene_name.mouse = "Gene") %>% left_join(mus.musculus.gene2ens)
maor_nof_cell_tdp_vs_ctrl.genelist = maor_nof_cell_tdp_vs_ctrl %>% select(gene_name.mouse, stat) %>% drop_na() %>% distinct(gene_name.mouse, .keep_all = TRUE) %>% tibble_to_matrix(stat, row_names = "gene_name.mouse")
# PROGENy
net_progeny.mouse <- get_progeny(organism = 'mouse', top = 100)
maor_nof_cell_tdp_vs_ctrl.progeny_scores = run_wmean(maor_nof_cell_tdp_vs_ctrl.genelist, net=net_progeny.mouse, .source='source', .target='target', .mor='weight', times = 1000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.003 ~ "****", p_value < 0.009 ~ "***",p_value < 0.01 ~ "**",p_value < 0.05 ~ "*", TRUE ~ ""))
maor_nof_cell_tdp_vs_ctrl_progeny.pathwayActivity <- ggplot(maor_nof_cell_tdp_vs_ctrl.progeny_scores, aes(x=reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred", mid = "whitesmoke", midpoint = 0) +
theme_oz(xaxis_arrow = FALSE) + theme(axis.text.x = element_text(angle = 90), legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = "TDP-43 over-expression", x = "", y = "TDP-43 overexpresison vs CTRL\nNormalised enrichment") +
stat_pvalue_manual(mutate(maor_nof_cell_tdp_vs_ctrl.progeny_scores, y.position = 9.5, group1=source, group2=source), label = "p.signif", xmin = "source", xmax = NULL, size = 4, hide.ns = TRUE)
maor_nof_cell_tdp_vs_ctrl_progeny.pathwayActivity
# Dorothea
net_dorothea.mouse <- get_dorothea(organism='mouse', levels=c('A', 'B', 'C'))
maor_nof_cell_tdp_vs_ctrl.dorothea_scores = run_wmean(maor_nof_cell_tdp_vs_ctrl.genelist, net=net_dorothea.mouse, .source='source', .target='target', .mor='mor', times = 1000, minsize = 5) %>% filter(statistic == 'norm_wmean') %>% mutate(p.signif = case_when(p_value < 0.05 ~ "*", TRUE ~ ""), sig = case_when(p_value < 0.135 & abs(score) < 4 ~ "sig", p_value < 0.135 & abs(score) >= 4 ~ "sig_strong", p_value >= 0.135 ~ "non_sig"), direction = ifelse(score > 0, "up", "down"), class = paste(sig, direction))
maor_nof_cell_tdp_vs_ctrl_dorothea.contrast_acts.volcano = ggplot(maor_nof_cell_tdp_vs_ctrl.dorothea_scores, aes(x = score, y = -log10(p_value))) + geom_point(aes(colour = class)) +
scale_colour_manual(values = c("non_sig up" = "gray", "non_sig down" = "gray", "sig up" = "#F36F6F", "sig_strong up" = "#EB4445", "sig down" = "#A6CEE3", "sig_strong down" = "#79B1D3")) +
labs(y = expression(-log[10]~P~value), x = "Normalised enrichment TDP-43 overexpression vs CTRL", title = "TDP-43 overexpression TFs") +
guides(colour = "none") +
theme_oz() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), panel.border = element_blank(), axis.ticks = element_line(colour = "black") ) +
geom_text_repel(fontface = "italic", data = filter(maor_nof_cell_tdp_vs_ctrl.dorothea_scores, source %in% "Tp53"), aes(label = source), max.overlaps = Inf, min.segment.length = unit(0, "lines"), size = 3) +
geom_vline(xintercept = 0, linetype = 'dotted') + ylim(c(0,2.8))
maor_nof_cell_tdp_vs_ctrl_dorothea.contrast_acts.volcano